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Abstract 

An algorithm for Monte Carlo simulations is proposed in which the param- 
eter controlling the strength of the transition becomes a dynamical variable 
and in which efficient transitions are achieved by cluster steps. Numerical 
results for the Potts model demonstrate the advantages of the method. 
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To overcome slowing down is crucial in computer simulations of phase transitions. 
For second-order phase transitions in a number of systems critical slowing down can 
be drastically reduced by cluster algorithms as introduced by Swendsen and Wang 
IQ. There the acceleration is achieved by the nonlocal update of larger structures. 

In the case of first-order transitions one encounters a different problem. On fi- 
nite lattices one has an exponentially fast suppression of the tunneling between 
metastable states of the system with increasing lattice size. To reduce this type of 
slowing down the multicanonical Monte Carlo algorithm has been proposed 0, ^ 
which uses enhancement of the suppressed configurations. The canonical distribu- 
tion is then obtained by an appropriate reweighting [Q. 

The application of cluster algorithms to first-order transitions leads only to little 
improvement ^. The obvious reason for this is that clusters can lead to some 
acceleration related to extended structures at not too large g, however, not with 
respect to the tunneling between order and disorder. 

To overcome slowing down in systems with a rough free-energy landscape the 
method of simulated tempering |^ has been proposed and applied to the random- 
field Ising model. In this method temperature values are related to a dynamical 
variable and the fact is utilized that at lower (3 the free-energy barriers are lower. 

In the present paper we propose to make the parameter which controls the strength 
of a first order transition a dynamical variable. This allows to proceed from a state 
at large strength to ones at lower strength where the transitions become easier. If in 
particular there is a range where the transition gets of second order, in this way the 
tunneling and the related suppression can be completely circumvented. It is crucial 
then to avoid also the critical slowing down related to extended structures. This is 
achieved by combining the procedure with cluster updates. 

We demonstrate in the case of the Potts model in two dimensions, using a 
square lattice with periodic boundary conditions, that this general scheme works 
very efficiently. The parameter which has to become dynamical in this case is q, the 
number of the degrees of freedom of the spins. 
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Figure 1 shows the distribution P{E, q) for energies and values q which is obtained 
for (3l (defined below). It is useful to illustrate the principle of our approach. Instead 
of tunneling across the valley we allow travelling along the peaks. 

Starting at an ordered state at large strength of the transition the cluster updates 
will mainly induce transitions between ordered states. Travelling down to lower 
strength they will increasingly break up larger structures. From low strength then 
an easy way is open to disorder at large strength. Further, also all the steps between 
different values of the strength variable will contribute to decorrelation. 

Before describing our update scheme it appears appropriate to formulate the 
general principle of making a parameter dynamical. A probability distribution 
fiq{a) of spin configurations {cr} which depends on a parameter q can be reinter- 
preted as the conditioned probability fl{q; a) for getting a configuration {cr} given 
a value of q. Then prescribing a distribution fl{q) one obtains the joint distribution 
cr) = fi{q)fi'{q', cr) for which the updates are to be performed. 

Our update scheme has three steps. The first one generates bond configurations 
{n} according to Swendsen and Wang and is described here by the conditioned 
probability A{a,q;n). The second step is a Metropolis update for the marginal 
distribution n) = X]{<t} /^('?) cr)^(<^) Q'! ^) realizing the conditioned probability 
T{n, q; q') for arriving at a new value q' . In the third step the new spins are ob- 
tained with the conditioned probability A{n,q']a') = ii{q' , a')A{a' ,q';n) / jl{q' ,n) 
which corresponds to the special choice of Swendsen and Wang where the new clus- 
ter spins are independent of the old ones. The generalized transition matrix thus 
becomes 

W{a,q-a\(^)=Y.A{a,q-n)T{n,q-q')A{n,q'-a') . (1) 

It leaves fi invariant if T leaves ft invariant, and it satisfies detailed balance with /j, 
if T does so with ft, which follows immediately from the definitions. 

In our update scheme cluster generation in addition to reducing critical slowing 
down also has the important virtue to allow technically simple steps between differ- 
ent g-values. The reason for this is that spins are not involved and the dependences 
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on the bond configurations are simple in jl{q,n). Using (i{q\ a) = fiq{cr) = Z 
with H = — ^(Ti,cTj and Z = J2{a} e"^"^, one gets its exphcit form 

j2{q,n) = /i(g)Z-i(g)(e^(«) - (2) 

where Nb is the number of bonds and Nc the number of clusters of the bond config- 
uration {n}. 

The simulations have been performed for two sets of g-values, q = A, ... ,7 and 
g = 2, . . . , 10 . Detailed balance for the Metropolis steps described by T(n, g; q') has 
been guaranteed using the symmetric proposal matrix Tp(g; q') = + + 

^g.toVte + Sg,qmaJq',qmaJ and the acceptance matrix mm{l,fi{n,q')/fi{n,q)). 

The distribution fl{q) has been choosen to be (approximately) constant within the 
sets of g- values considered. In this context also the g-dependence of Z^^{q) in @ 
is to be known. By requiring the number of sweeps to be the same for each of the 
g- values in the respective set it has been determined numerically in an iterative way. 
Table 1 shows the relative values of Z~^{q)Ji{q) obtained for the set g = 4, . . . , 7 
requiring constancy of /2(g) better than 4% . 

For the values of /3(g) two choices have been considered. The first one, denoted 
by /5i(g) here corresponds to equal height of the two peaks in the distribution of 
energies for g > 4 or to the maximum of the specific heat which is used for g < 4. 
Table 2 gives the values of /^^(g) used for the set g = 4, . . . , 7. The second /5-value 
considered is the infinite- volume critical value (3c =ln(y/g + 1). 

The observables measured for all g are the energy E and the order parameter 
M = {q max(a)(Ai'a) — l)/(g — 1) with Na = V"^ J^i^aia ■ For both of them in each 
case also exponential and integrated autocorrelation times have been determined. 
Further, for comparison also simulations for other algorithms have been performed 
and in addition tunneling times determined. The statistics collected has been in the 
range of 3 x 10^ to 2 x 10^ sweeps for each g and all algorithms investigated. 

There are results on other algorithms in literature for g = 7 |10|, and for 
g = 10 1^, P). We performed simulations for the heat bath algorithm (H), the 



method of Swendsen and Wang (SW), and the multicanonical heath-bath algorithm 
(MH) ourselves to get autocorrelation times for a precise comparison (rather than to 
have to rely on tunneling times only available in most cases). Because the Metropolis 
algorithm has been found to be rather slow it has not been considered here. SW 
turns out to be slightly superior to H if one accounts for CPU times, too. 

Tunneling times |l^ are only useful at larger q where they get dominant in the 
autocorrelation times. We observe that they depend on the particular definition and 
differ from the autocorrelation times, which for MH and DQ persists even at larger 
L. 

Table 3 gives the average stay times for each q (the average number of sweeps 
spent at a particular q before leaving it) obtained for [3l and the set g = 4, . . . , 7 . 
It is seen that the travelling goes relatively fast. For f3c and for the set g = 2, . . . , 10 
the stay times are very similar. 

In Table 4 the exponential autocorrelation times for (3l and the set g = 4, . . . , 7 are 
shown. The errors are statistical ones. In the cases of (3c and of the set g = 2, . . . , 10 
the results are very similar. The results obtained for the order parameter within 
errors agree well with those for E (listed in the Table). Thus the measurement of 
both observables provides a valuable check. 

Figures 2 and 3 give the exponential autocorrelation times for g = 7 and g = 10, 
respectively, which we obtained for the heat bath algorithm (H), the method of 
Swendsen and Wang (SW), the multicanonical heath-bath algorithm (MH), and our 
approach with dynamical g (DQ) for the sets g = 4, . . . , 7 and g = 2, . . . , 10 . The 
time scale is in sweeps of the respective algorithm in each case. 

For a comparison firstly the CPU times related to these scales are to be considered. 
Because the fraction of time for the Metropolis step is small (of the order of 4%) as 
compared to that of the cluster steps in DQ, these scales for SW and DQ are about 
the same. For H and MH the underlying sweeps are about a factor 2 slower. 

To compare secondly one has to discuss that in DQ time is also spent at the other 
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values of q. With the extreme view that only one of the g-values is of interest - and 
all other results obtained at no extra cost are ignored - one could multiply by 4 and 
by 9 in the cases of the sets q = 4, ... ,7 and g = 2, . . . , 10 , respectively (noting 
that according to the choice of fi{q) equal times are spent at different q). Even then 
one gets an improvement as compared to MH of up to factors 7 and 3 for g = 7 and 
q = 10, respectively, in the range considered With the more realistic view that 



all results are of interest the gain is more than an order of magnitude. 



For g = 7 by MH as compared to H there is only a modest gain |]T0| of about a 



factor of 2. Recently by a multicanonical demon algorithm |1T| some improvement 
upon MH has been achieved on larger lattices. For DQ on the set g = 4, . . . , 7 we 
find an improvement with respect to both of these multicanonical algorithms which 
even after a multiplication by 4 as mentioned above is still about a factor 6 . Thus 
by DQ in any case one improves more than an order of magnitude as compared to 
the conventional case. 

For g = 10 MH gets more efficient, gaining more than an order of magnitude in 
the range considered here as compared to H. However, DQ remains still superior 
even after the multiplication by 9 mentioned above. In addition, the multiplication 
factor can be reduced by using the set g = 4, . . . , 10 instead of the set g = 2, . . . , 10. 
This is suggested by the stay times indicating that it is not worthwhile to go down 
to g = 2. This feature is also apparent from the comparison of the results for DQ 
for g = 7 obtained on the sets g = 4, . . . , 7 and g = 2, . . . , 10 where the advantage 
of not going down amounts to a factor of about 2. 

In this context it should be noted that there are a number of possibilities for 
further improvements. Some have already been briefly tested with positive results. 
For example, instead of g-steps by one one can use steps by two. To accelerate 
certain transitions one can make f3 dynamical, too. Instead of keeping /i(g) constant 
it may be adjusted to get more favorable stay times. Modification of the Metropolis 
step offers various optimizations. 

For g = 4 and 3, where the phase transition gets of second order, comparing with 



6 



integrated autocorrelation times obtained for SW [O one observes still reductions 



by factors of about 4 and 3, respectively. For q = 2 comparing with exponential 
autocorrelation times for SW one gets a corresponding factor of about 1.4 . 
To these reductions the additional tunneling possible at not too large L contributes 
which does not persist at very large L. 

From Figures 2 and 3 it is seen that for DQ there is some small exponential- 
like contribution present in the autocorrelation times. The obvious reason for this is 
that at not too large L there are still sizable additional contributions from tunneling 
which die out at larger L. At very large L only travelling along the peaks is possible 
and the true asymptotic behavior occurs. 

At the upper end of the range considered for DQ the increase of autocorrelation 
times is smaller than for MH which gives first hints about the ultimate asymptotic 
behavior. However, to conclude about the ultimate behavior much larger lattices 
are needed and it is to be noted that then details as, for example, the appropriate 
adjustment of fl{q) become important. From the practical point of view on the 
lattices which can be reached there remains in any case substantial gain for DQ. 

The proposed concept has been realized here with the simplest choices of details 
in order to check its working in practice. Clearly the next task is to optimize 
these choices. In any case the present results demonstrate that the concept works 
surprisingly well such that it deserves further development and also application to 
other systems. 
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Table 1 

Values of {Z-\q)fi{q)) / {Z-\7)ji{7)) . 



L 


g = 6 


g = 5 


g = 4 


12 


7.23 X 10"^ 


5.43 X 10-1 


5.346 X 10-1 


16 


4.748 X 10-1 


2.629 X 10-1 


1.925 X 10-1 


24 


1.726 X 10-1 


3.669 X 10-2 


1.677 X 10-2 


34 


2.857 X 10-2 


1.18 X 10-^ 


1.95 X 10-^ 


50 


5.17 X 10-^ 


5.38 X 10-^ 


1.52 X 10-^ 



Table 2 



Numerical values for (3l ■ 



L 


q = 7 


q = 6 


q = 5 


q = 4 


12 


1.2725 


1.2158 


1.1495 


1.0708 


16 


1.2806 


1.2238 


1.1582 


1.0792 


24 


1.2872 


1.2309 


1.1656 


1.0879 


34 


1.29005 


1.23407 


1.16915 


1.0918 


50 


1.29178 


1.23607 


1.17145 


1.0947 
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Table 3 

Average stay times for Pl (with errors smaller than 1%). 



L 


q = 7 


q = 6 


q = 5 


q = 4 


12 


2.48 


1.27 


1.32 


2.70 


16 


2.70 


1.39 


1.47 


3.08 


24 


3.22 


1.69 


1.88 


4.02 


34 


4.19 


2.24 


2.66 


5.85 


50 


6.84 


3.81 


5.04 


12.8 



Table 4 



Exponential autocorrelation times ioi (3l . 



L 


q = 7 


q — 6 


g = 5 


g = 4 


12 


15.1(2) 


11.5(1) 


9.04(6) 


8.18(6) 


16 


21.2(4) 


15.8(2) 


12.5(1) 


10.9(1) 


24 


36.8(6) 


26.7(3) 


20.4(2) 


16.6(2) 


34 


64.0(6) 


44.7(5) 


31.8(3) 


24.1(3) 


50 


126(2) 


79(1) 


50(1) 


37.0(4) 
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Figure Captions 

FIG. 1. Distribution P{E, q) for (shown for L = 34). 

FIG. 2. Comparison of exponential autocorrelation times for q = 7 and 

the algorithms H, MH, DQ on set g = 4, . . . , 7, and DQ on set 

g = 2,...,10. 

FIG. 3. Comparison of exponential autocorrelation times for q = 10 and 
the algorithms SW, H, MH, and DQ on set g = 2, . . . , 10. 
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